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We present a continuum theory to describe elastically induced phase transitions between coherent 
solid phases. In the limit of vanishing elastic constants in one of the phases, the model can be used to 
describe fracture on the basis of the late stage of the Asaro-Tiller-Grinfeld instability. Starting from 
a sharp interface formulation we derive the elastic equations and the dissipative interface kinetics. 
We develop a phase field model to simulate these processes numerically; in the sharp interface limit, 
it reproduces the desired equations of motion and boundary conditions. We perforin large scale 
simulations of fracture processes to eliminate finite-size effects and compare the results to a recently 
developed sharp interface method. Details of the numerical simulations are explained, and the 
generalization to multiphase simulations is presented. 

PACS numbers: 62.20.Mk, 46.50.+a, 81.40.Np 
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I. INTRODUCTION 

The propagation of cracks is very important for many 
applications and a central topic in physics and materials 
science. The most fundamental basis of understanding 
fracture traces back to Griffith [l[ ; according to his find- 
ings, the growth of cracks is determined by a competition 
of a release of elastic energy and a simultaneous increase 
of the surface energy if a crack extends. Although much 
progress has been made in understanding the striking fea- 
tures of cracks [2[, the mechanisms which determine the 
dynamics of crack propagation are still under heavy de- 
bate. A typical description of cracks starts on the atomic 
level and interprets the propagation by successive break- 
ing of bonds; it is obvious that the theoretical predictions 
significantly depend on the underlying empirical models 
of the atomic properties (see for example [3|). A rather 
complementary approach takes into account effects like 
plasticity, which can lead to extended crack tips (finite tip 
radius ro) Recent experimental investigations of frac- 
ture in brittle gels 0] possibly reveal macroscopic scales. 
It is obvious that under these circumstances a full mod- 
eling of cracks should not only determine the crack speed 
but also the entire crack shape and scale self-consistently. 

During the past years, phase field modeling has 
emerged as a promising approach to analyze fracture 
by continuum methods. Recent phase field models go 
beyond the microscopic limit of discrete models, and 
encompass much of the expected behavior of cracks 
[1,0,1,1; However, a significant feature of these descrip- 
tions is that the scale of the growing patterns is always set 
by the phase field interface width, which is a purely nu- 
merical parameter and not directly connected to physical 
properties; therefore these models do not possess a valid 
sharp interface limit. Alternative descriptions, which are 
intended to investigate the influence of elastic stresses on 
the morphological deformation of surfaces due to phase 
transition processes, are based on macroscopic equations 



of motion. But they suffer from inherent finite time sin- 
gularities which do not allow steady state crack growth 
unless the tip radius is again limited by the phase field 
interface width 10 . Very different approaches which are 
not based on a phase field as order parameter introduce 
a tip scale selection by the introduction of complicated 
nonlinear terms in the elastic energy for high strains in 
the tip region [ll|, requiring additional parameters. 

Recently, we developed a minimum theory of fracture 
[l2| which is only based on well-established thermody- 
namical concepts. This is also motivated by experimen- 
tal results showing that many features of crack growth 
are rather generic [l3j ; among them is the saturation of 
the steady state velocity appreciably below the Rayleigh 
speed and a tip splitting for high applied tension. This 
theory describes crack growth as a consequence of the 
Asaro-Tiller-Grinfeld (ATG) instability [IJ] in the frame- 
work of a continuum theory. Mass transport at the ex- 
tended crack tips can be either due to surface diffusion 
or a phase transformation process. The latter has been 
investigated numerically by pha se field simulations [l5j 
and sharp interface methods [16|. It turned out that the 
phase field simulations were still significantly influenced 
by finite size effects and insufficient separation of the ap- 
pearing length scales, and therefore the results did not 
coincide. One central aim of the current paper is there- 
fore to carefully extrapolate new phase field results ob- 
tained by large-scale computations. As we will show, we 
then get a very convincing agreement of the approaches. 
Also, we explain details of the phase field method and 
the underlying sharp interface equations in more detail. 

The paper is organized as follows: First, in Section 
HI! we introduce the sharp interface equations to de- 
scribe crack propagation as a phase transformation pro- 
cess. The basic selection mechanisms for crack growth 
are reviewed in Section IIII1 We introduce a phase field 
description to solve the arising moving boundary problem 
in Section [TV] We demonstrate the numerical separation 
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of lengthscales and obtain results which are in excellent 
agreement with sharp-interface predictions (Section fVf . 
This part contains important new results concerning the 
underlying continuum theory of fracture. In Section I VII 
we briefly explain how the model can be extended to sys- 
tems consisting of multiple phases. Detailed derivations 
of the sharp interface equations are given in Appendix 
lAl Since coherency of two solid phases leads to an unex- 
pected expression for the chemical potential, its relevance 
is analyzed specifically for the motion of a planar inter- 
face (Appendix [Bj. In Appendix [Cl it is demonstrated 
that the presented phase field model recovers this effect 
in the sharp interface limit. Finally, details of the numer- 
ical implementation of the phase field model are given in 
Appendix [Dl 



II. SHARP INTERFACE DESCRIPTION 

The fracture process in [l5| is interpreted as a first 
order phase transition from the solid to a "dense gas 
phase" , driven by elastic effects. More generally, we 
investigate the transition between two different solid 
phases. In the limiting case that one phase is infinitely 
soft, crack propagation can be studied. A central simpli- 
fication is due to the assumption of equal mass density p 
in both phases and the condition of coherency at the in- 
terface, i.e. the displacement field m is continuous across 
the phase boundary, 
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where the upper indices indicate the different phases. 
The strain is related to the displacement field Ui by 
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Strain and stress Oik are connected through Hooke's law; 
for the specific case of isotropic materials, it reads 
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for each phase. Here, E and v are elastic modulus and 
Poisson ratio, respectively. We note that eigenstrain con- 
tributions due to different unit cells of the phases are not 
considered here for brevity, but they can easily be intro- 
duced [13] ■ For simplicity, we assume a two-dimensional 
plane-strain situation. 

All following relations are obtained in a consistent way 
from variational principles, and this is described in detail 
in Appendix[X] The equations of dynamical elasticity are 
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for each phase. On the interface, we obtain the expected 
continuity of normal and shear stresses 
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FIG. 1: Geometry of the phase transition scenario. Phase 
transitions between phases 1 and 2 are possible and lead to 
interface motion with local normal velocity v n . The volumes 
of the two phases and the interface A(t) are therefore time- 
dependent. 



Here, the index n denotes the normal direction of the 
interface, with the perpendicular tangential direction t 
(see Fig. [I}. We mention that this equation holds only 
for the specific conditions of equal mass densities and 
coherency at the interface. In other cases one obtains 
more general relations for momentum conservation at the 
interface, which also involve the interface normal velocity 

The elastic contribution to the chemical potential at 
the interface for each phase is 
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Since the chemical potential has the dimension energy 
per particle, we introduced the atomic volume Cl. It is 
quite remarkable that the normal and shear contributions 
enter into the expression with negative sign, in contrast 
to the natural expectation p e i — d<Jik€ik/2, which is the 
potential energy density. The reason for this modification 
is the coherency constraint which has to be fulfilled at 
the interface. An illustrative example to understand this 
unexpected expression for the chemical potential is given 
in Appendix [B] However, this effect is only important 
for solid-solid transformations. For crack propagation, 
where we assume that the new phase inside the crack is 
infinitely soft, the normal and shear stresses vanish at 
the interface, and therefore the discrepancy between the 
chemical potential ^ and the naive guess p e i disappears. 

We would also like to mention that the equality of the 
mass density and the coherency leads to the absence of ki- 
netic energy contributions to the chemical potential. The 
reason is that such a contribution is continuous across 
the interface; finally, only the chemical potential differ- 
ence fj,^ — p^ enters into the equation of motion for the 
interface, and therefore such a term would cancel. Nev- 
ertheless, we note that the kinetic energy contribution 
would enter into the chemical potential with negative 
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sign, i.e. —piif, because in the Lagrangian the kinetic 
and potential energy contributions appear with opposite 
sign. Kinetic contributions may play a role if instead of 
a phase transformation process, surface diffusion along a 
free boundary drives the evolution. 

For the motion of the interface, surface energy is also 
taken into account. Since it does not couple to the elastic 
terms, it simply gives an additional contribution to the 
chemical potential, 

with the local interface curvature k and the surface en- 
ergy density 7. The curvature is positive if phase 1 is 
convex. Then the motion of the interface due to a phase 
transition process is described by the local normal veloc- 
ity (D is a kinetic coefficient with dimension [D] = m 2 /s) 
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which is positive if phase 1 grows. The set of equations 
JI])-© describes the dynamics of the system. We point 
out that it leads to a complicated free boundary problem, 
and the arising interfacial patterns are self-consistently 
selected during this nonequilibrium process if external 
forces are applied to the system. For an initially almost 
flat interface between a soft and a hard phase, which 
is nonhydrostatically strained, first the ATG instability 
develops: Long wave morphological perturbations lead to 
a decrease of the total energy and the formation of deep 
notches, similar to cracks (see e.g. [l0|, EH)- As we have 
shown in [l2|, EH EH , it is essential to include the inertial 
contributions, because otherwise a steady state growth of 
these cracks is impossible, and the system collapses into 
the finite time cusp singularity of the ATG instability. 



III. CRACK PROPAGATION: SELECTION 
PRINCIPLES 

In the case that one phase is infinitely soft, crack prop- 
agation can be studied. Here, growth of the crack is based 
on a phase transition of the solid matrix to a "dense gas 
phase" which has the same density as the solid. In this 
sense, it is similar to other models of fracture based on a 
non-conserved order parameter [1, @, [1] • The crucial dif- 
ference is that the current model is based on well-defined 
sharp interface equations, and therefore the predictions 
do not depend on inherently numerical parameters like 
a phase field interface width. However, numerically, it 
requires a tedious separation of scales to obtain these 
results; this is described in the next section. 

Understanding fracture as a phase transition process 
offers many numerical advantages, as phase field models 
can be derived to solve this moving-boundary problem. 
We point out that the underlying selection principles 
which allow a steady state crack growth with propaga- 
tion velocities well below the Rayleigh speed, tip blunting 



and branching for high driving forces are rather generic 
and are similarl y v alid for models with conserved order 
parameters. In {12}, we derived the similar equations of 
motion if instead of a phase transition, mass transport is 
due to surface diffusion along the free crack boundary. 

In the latter case, the elastic boundary conditions are 
replaced by 
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and the chemical potential becomes 
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It differs from the expression ([6]) first by the elastic en- 
ergy density, because no coherency constraints have to be 
fulfilled here. Second, the kinetic energy density appears 
here, because it does not cancel in the derivation from 
the Lagrangian. The equation of motion for the interface 
is replaced by 
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with the surface diffusion coefficient D^ sd \ 

In both cases, stresses on the boundary of the crack 
tip with finite radius r scale as 



Kr, 
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and the curvature behaves as k ~ l/^o- Therefore, 
all contributions to the chemical potentials scale like 

— 1/2 

fi ~ r Q ' , and this is ultimately the reason for the cusp 
singularity of the Grinfeld instability and the impossibil- 
ity of a steady-state crack growth, if only static elasticity 
is taken into account: Then, the equations of motion ([7]) 
or (fT0|) can be rescaled to an arbitrary tip radius which is 
not selected by the dynamical process. The explanation 
is that the linear theory of elasticity and surface energy 
define only one lengthscale, the Griffith length, which 
is macroscopic, but do not provide a microscopic scale 
which allows the selection of a tip scale. Formally, the 
equations of motion depend only on the dimensionless 
combinations vro/D for the phase transition dynamics 
and 

vr 3/ D (sd) for 

surface diffusion; the radius ro and the 
steady state velocity v therefore cannot be selected sep- 
arately; any rescaling which maintains the value of the 
product would therefore describe another solution. The 
situation changes if inertial effects are taken into account, 
which is reasonable for fast crack propagation. Then ad- 
ditionally the ratio v/vr (vr is the Rayleigh speed) ap- 
pears in the equations of motion, and therefore a rescal- 
ing is no longer possible. Instead, D/vr for the phase 
transition dynamics and (£)( sd ) /vr) 1 ^ 3 for surface diffu- 
sion set the tip scale. Thus, we conclude that fast steady 
state growth of cracks is possible if inertial effects are 
taken into account. More formal analyses also including 
rigorous selection mechanisms due to the suppression of 
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growing crack openings far behind the tip are given in 
|15| for phase transition processes and in [12] for surface 
diffusion. 

This analysis shows that the selection principles which 
allow a fast steady state growth of cracks are similar 
for the simple phase transition process studied here and 
surface diffusion. The latter does not require the intro- 
duction of a dense gas phase inside the crack and obeys 
conservation of the solid mass itself. Even more, surface 
diffusion can be understood in a generalized sense as plas- 
tic flow in a thin region around the extended tip which 
can be effectively described in the spirit of a lubrication 
approximation. Therefore, many general statements ob- 
tained for the phase transition dynamics can also be used 
for crack growth propelled by surface diffusion. The lat- 
ter is more tedious to implement numerically, since the 
equation of motion (fT0|) is of higher order [19[ . 

For both mechanisms, a tip splitting is possible for 
high applied tensions due to a secondary ATG instabil- 

— 1/2 

ity: Since a ~ Kr in the tip region and the local 
ATG length is Lq ~ Ej/a 2 , an instability can occur, 
provided that the tip radius becomes of the order of the 
ATG length. In dimensionless units, this leads to the 
prediction A spUt ~ 0(1). 



IV. PHASE FIELD MODEL 

To describe systems with moving boundaries accord- 
ing to the equations of motion developed above, we im- 
plemented a phase field model. Let <f> denote the phase 
field with values <f> = 1 for phase 1 and <f> — for phase 
2. The energy density contributions are 



f el =ii{cj>)el +A(0)(e i2 ) 2 /2 
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for the elastic energy, with the interpolated shear modu- 
lus and Lame coefficient 

K<P) = h{4>)^ + (\-h{4>))^ 2 \ (13) 

\{<f>) = h(<f>)\W + (1 - h(cf>))\W (14) 



where 
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interpolates between the phases, and the superscripts de- 
note the bulk values. The surface energy is 
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with the interface width £. Finally, 
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is the well-known double well potential. Thus the total 
free energy is given by 
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The elastodynamic equations are derived from the free 
energy by variation with respect to the displacements m, 
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and the dissipative phase fields dynamics follows from 
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F= dV (f el + f s + f dw ) . 



It has been shown in [1C| that in the quasistatic case, the 
above equations lead to the sharp interface equations Q- 
([7]) if the interface width £ is significantly smaller than 
all physical lengthscales present in the system. In Ap- 
pendix [Cl it is illustrated that this model also correctly 
incorporates the modification of the chemical potential 
© due to the coherency constraint. Details of the nu- 
merical implementation are given in Appendix [D] 



V. PHASE FIELD MODELING OF CRACK 
PROPAGATION 

The central prediction of this theory of fracture is that 
a well-defined steady-state growth with finite tip radius 
and velocities appreciably below the Rayleigh speed is 
possible. This also cures the problem of the finite-time 
cusp singularity of the Grinfeld instability. These pre- 
dictions have been confirmed by phase field simulations 
[TH and sharp interface methods [l6[ which are based on 
a multipole expansion of the elastodynamic fields. Sur- 
prisingly, it turned out that the obtained results seem 
to differ significantly: In particular, the sharp interface 
method predicts a range of driving forces inside which 
the velocity of the crack is a monotonically decreasing 
function. Here, we demonstrate that the discrepancy of 
results is due to finite size effects of the previous phase 
field results [TB|, and that by careful extrapolation of 
large-scale simulations, a coinciding behavior is obtained. 

We investigate crack growth in a strip geometry with 
fixed displacements at the upper and lower grip. The 
multipole expansion technique [l6j is designed to model 
a perfect separation of the crack tip scale D/vr to the 
strip width L: In most real cases, crack tips are very 
tiny, and therefore it is theoretically desirable to describe 
this limit. For the phase field method, however, a finite 
strip width L is necessary, and a good separation of the 
scales therefore requires time-consuming large-scale cal- 
culations. We typically use strip lengths 2L and shift 
the system such that the tip remains in the horizontal 
center. This allows to study the propagation for long 
times until the crack reaches a steady state situation. 
Apart from this finite size restriction, we had to intro- 
duce the interface width £ as a numerical parameter, and 
the phase field method delivers quantitative results only 
in the limit that all physical scales are much larger than 
this lengthscale. The latter has to be noticeably larger 
than the numerical lattice parameter Ax, but the results 
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FIG. 2: Crack shapes for different scale separations D/vr^ 
and fixed ratio Lvr/D — 1 1 .03; the aspect ratio of the system 
is 2 : 1. The driving force is A = 1.4. By improvement of 
the separation, the crack opening is reduced, and finally the 
boundaries become straight parallel lines. 




FIG. 3: First step of the extrapolation procedure for the di- 
mensionless velocity v/vr. The system size L/£ is increased 
and the ratio D/^vr is kept fixed for each curve. For each 
ratio, an extrapolated velocity v^/vr corresponding to an in- 
finite system size is obtained, as indicated by the dashed lines. 
The driving force is A = 1.4. 



show that the choice £ = 5 Ax is sufficient. We therefore 
have to satisfy the hierarchy relation 
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which is numerically hard to achieve. We developed a 
parallel version of the phase field code which is run on 
up to 2048 processors, with system sizes up to 8192 x 
4096 ■ (Ax) 2 . All computations are performed on the su- 
percomputers JUMP and JUBL operated at the Research 
Center Jiilich. 

For the strip geometry, a dimensionless driving force 
A is defined as 



S 2 (X + 2^) 
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with S being a fixed displacement by which the strip is 
elongated vertically. The elastic constants of the new 
phase inside the crack are zero. The value A = 1 cor- 
responds to the Griffith point. All calculations are done 
with Poisson ratio v = 1/3. 

In [l6|], it was shown that close to the Griffith point, 
dissipation free solution exists in the framework of the 
model: In this regime 1 < A < 1.14 an additional micro- 
scopic length scale is needed to select the small tip radius 
which is no longer determined by the ratio D/vr. This 
can here be mimicked by the phase field interface width 
and was already done in flEl ]. Here, we focus on the more 
interesting regime of higher driving forces, but still below 
the threshold of instability. Typical crack shapes in the 
vicinity of the tip are shown in Fig. [21 

To fulfill the scale separation (12~1]) , we perform a dou- 
ble extrapolation of the obtained steady state velocities 
vl,£ (the subscripts indicate the additional non-resolved 



length scale dependencies). In the first step, we ex- 
tend the simulations to an infinite system size. There- 
fore, we decrease the ratio £/£ — > for fixed tip scale 
ratio D/£vr. This step is demonstrated in Fig. [3] for 
A = 1.4. Here, the dimensionless propagation velocity 
Vl£ = vl^/vr is plotted as function of the inverse square 
root of the system size (^/L) 1 / 2 . In this representation, 
the data for the larger systems can be extrapolated lin- 
early to infinite system sizes, since we numerically get a 
scaling 
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for large systems, £/L <C 1, with a constant a > for 
each curve. Since the separation of D/vr to £ is still 
imperfect, the extrapolated values «^(A, D/vr£) do not 
yet cumulate to a single point, and therefore a second 
extrapolation step is necessary. 

Hence, in Fig.[U the dependence of the velocity v^/vr 
on the separation parameter vr^/D for A = 1.4 is shown. 
The extrapolated values from Fig. [3] are used, and we 
obtain a scaling 



*(A,iL)= e (A)-^ 



^ (24) 



with a constant /3 > and the dimensionless sharp inter- 
face limit velocity v — v/vr. 

This tedious procedure was performed for several driv- 
ing forces, and in Fig. [5] the comparison to the multi- 
pole expansion method 16] is shown. The agreement 
of the results which are obtained from completely differ- 
ent methods is very convincing. The small deviation for 
A = 1.8 is due to the fact that this value is already close 
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FIG. 4: Second extrapolation step to obtain the sharp in- 
terface velocity v. The extrapolated velocities obtained from 
Fig. [3] are plotted as function of the scale separation param- 
eter va/^D. In this example A = 1.4 is used. 
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FIG. 5: Comparison of the steady-state crack velocity ob- 
tained from the multipole expansion technique [l6j] and the 
extrapolated value from phase field simulations. 



to the threshold of the tip-splitting instability which can- 
not be captured by the multipole expansion method. In 
particular, we find evidence for the prediction that the 
steady state velocity decays weakly with increasing driv- 
ing force, which might be an artefact of the model. 

For higher driving forces, we observe tip-splitting in 
the phase-field simulations (see Fig. [S]). The onset of 
this irregular behavior depends sensitively on the system 
size, because in relatively small systems, the branches 
of the crack cannot separate since they are repelled by 
the boundaries. Therefore, the steady state growth is al- 
ways stabilized by finite size effects. On the other hand, 
initial conditions can trigger an instability, and then a 
long transient is required to get back to steady state so- 
lutions. Despite these restrictions, we are still able to 
make the prediction that the threshold of splitting obeys 
A S put < 1.9 in the phase field model. It is in agreement 



FIG. 6: Irregular tip splitting scenario for A = 1.9. We used 
Lvr/D = 44.2 and D/vr£, = 9.3; the aspect ratio of the 
system is 2:1. Time is given in units D/v R . The thickness 
of the interface is the phase field interface width, indicating a 
good separation of the scales. 



with the conjecture that branching occurs as soon as the 
steady state tip curvature becomes negative, leading to 
the prediction A sp ii t « 1.8 [l6j |. 

The numerical determination of a characteristic crack 
width scale in the sharp interface limit is more difficult, 
and therefore we refrain from performing a double ex- 
trapolation procedure. The explanation is that if the 
soft phase inside the crack still possesses small nonva- 
nishing elastic constants, the equilibrium situation far 
behind the crack tip corresponds to a full opening of the 
crack, instead of the opening being of the order D/vr: 
As it is shown in Appendix [Bj the elastic energy is min- 
imized if the hard phase completely disappears. Small 
remaining elastic constants can be due to an insufficient 
separation of the scales D/vr and £, since according to 
Eqs. (H2J) and (fT4"|) . the elastic constants decay only ex- 
ponentially inside the crack, even if this soft phase has 
nominally vanishing clastic coefficients. Therefore, the 
crack opening is a weakly growing function of the dis- 
tance from the crack tip, and this slope becomes smaller 
with better scale separation, see Fig. [5] We point out that 
this opening is solely due to the phase transition process, 
and the shapes are drawn without elastic displacements 
which should be added to obtain the real shape under 
load. For example, the vertical displacement obeys the 
usual scaling u y ~ y\x\ for large distances \x\ from the 
tip. 

The same effect can be seen if we investigate solid-solid 
transformations towards a soft phase with small elastic 
constants. The Poisson ratios in both the surround- 
ing solid and the new inner phase are chosen equally, 
v = 1/3, but the bulk moduli differ by many orders of 
magnitude. The softer the inner phase becomes, the less 
the opening of the "crack" grows with increasing distance 
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FIG. 7: Solid-solid transformation in a strip. A very soft 
phase grows at the expense of a harder phase. Parameters 
are Lvr/D = 11.03 (vr is the Rayleigh speed of the hard 
phase), D/£vr = 9.27, the aspect ratio is 2 : 1. The Poisson 
ratio v = 1/3 is equal in both phases, and the driving force is 
A = 1.4. 



from the tip, see Fig. [71 Only very far away, the new 
phase fills the whole channel. 



VI. MULTIPHASE MODELING 

A simple approach to derive multiphase equations to 
describe systems consisting of more than two phases 
starts again from variational principles. The volume 
fraction of each phase is described by a field variable 
<pk, k = 1, . . . , N with N being the number of phases, 
<f> = (0i, . . . , 4>n). In the sharp interface limit, one phase 
field variable has the value one inside the bulk phases, 
and the others are zero. Their temporal evolution is given 
by 
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where we introduce a Lagrange multiplier to maintain the 
phase conservation, Y2k=i 4> k ~ 1- We redefine the diffu- 
sion coefficient D/j — > D, which is more appropriate to 
generalize to arbitrary interfacial energy coefficients jik 
between phases i and k. The additional factor 2 in the 
equation of motion above is chosen to recover the previ- 
ous phase field model in the case N — 2. The expression 
for the Lagrange multiplier is given by 
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The free energy F 
contributions: 



F e i + F s + F dw has the following 



F e j = 
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F s = fJ2 7« / - &V&) a dV, (28) 
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Here, the interpolated elastic constants are 

N N 

n(*) = YM)n®, A($)=y>(fc)AW (30) 
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with //W an d \W being the elastic constants of the indi- 
vidual bulk phases. Also, we have mutual interfacial en- 
ergy coefficients 7^ = 7^ for phase boundaries between 
i and j. Notice that third-phase contributions appear in- 
side two-phase boundaries which lead to a renormaliza- 
tion of the bare interfacial energies. This can be avoided 
by the addition of higher order terms to the multiwell- 
potential (f2"9")l . See Refs. [13, [SJ for a discussion of this 
issue. 

Variation gives explicitly: 



SF 



el 



0(pk 

86k 



h'{6 k ) (^ k el + X k {e u ) 2 /2) , 

N 
i=l 



SF, 



(I if 



12 " 

7- 2^ ™ 

> 1=1 



Similarly, the elastic equation of motion is just the gen- 
eralization of Eq. (fTi?]) 



pui 



dx k 



with 



CT ifc ($) = 2//(4>)e 4fc + \($)euS ik . 



(31) 



(32) 



This description is a direct generalization of Ref . [2pJ , and 
therefore it is known that it leads to appropriate sharp 
interface equations similar to Eqs. JT])-((7|) and contact 
angles as predicted by Young's law [21j ]. 



APPENDIX A: DERIVATION OF THE SHARP 
INTERFACE EQUATIONS 

Here we explain in detail how the equations of motion, 
the appropriate boundary conditions and the chemical 
potential, which is responsible for interface motion, are 
derived in a unique way from variational principles. We 
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assume that the two phases are coherent, i.e. the dis- 
placement field is continuous across the interface, and 
the mass densities are equal in both phases. Since we 
do not consider lattice strains or surface tension here, 
all elastic stresses arise from external forces. For sim- 
plicity, we assume a two-dimensional plane-strain situa- 
tion. Contributions due to surface energy do not couple 
to the elastic fields. Ultimately, they give only an addi- 
tive curvature-dependent term to the chemical potential, 
as incorporated in Eq. ([7]). We note that this expression 
was already obtained in [22| . but here we take also in- 
ertial effects into account, since the velocity of cracks is 
typically of the order of the sound speeds. 
The kinetic energy density is in both phases 



rr 1-2 



and the potential energy density reads 

r/(«) - !„-(")>) 

u ~ 2 *k t ik ■ 



(Al) 



(A2) 



Notice that certain components of the stress and strain 
tensors are in general discontinuous at the interface, as 
will be elaborated below. 

We assume the total volume V of the entire system 
to be constant in time and to be decomposed into two 
subvolumes V^(t) and V^ 2 '(t) of different solids. Up- 
per indices discriminate between the phases (see Fig. [TJ . 
The common interface A(t) := dV {1) {t) n <9U (2) (i) with 
normal n and tangential r is moving in time due to 
phase transitions, and consequently, the phase volumes 
are time-dependent as well. However, we do not yet spec- 
ify a concrete dynamical process here. 

The Lagrangian is defined as 

L(t) = J TdV - J U [1) dV - J U {2) dV, (A3) 

V Vi(t) V 2 (t) 

and the action is 

ti 

(A4) 



S = J L(t)dt, 



with arbitrary beginning and end times to and t\. 

We obtain the usual elastic equations by varying the 
action (|A4[) with respect to the displacement field for 
fixed interface positions. Thus we get 



ti 

SS = J dt 

to 



pUiSiiidV - 

V Vi(t) 



o^S^dV 



(2)c (2) „ r 



V 2 {t) 



dt 



f da (1) 

dV+ / — i^SuidV 

J ox k 

Vi(t) 



a^l^Suidr + 



(2) 



dxi- 



■SuidV 



A(t) 



V 2 (t) 



a^) 2) Suidr 



A(t) 



The first integral is integrated by parts, assuming as 
usual that the variations Sui vanish for to and t\. Since 
also the normal vectors of both phases are antiparallcl, 



n := n« 



-n (2) , thus 



CT;„(2) 



5S = 



dt 




-cr m , we get 



pui Su.idV 



■Vi(t) 




- pui ) Su.dV 
dx k ) 



Demanding vanishing variation 5S gives in the bulk the 
usual equations of motion 



0, 



dxk 



PUi, 



(A5) 



and on the interface we obtain the continuity of normal 
and shear stresses 



(A6) 



The next step is to calculate the change of the total 
energy when the interface moves in the course of time. 
This is done in three steps: First, we calculate the change 
of energy due to the time evolution of the elastic fields for 
fixed interface position. Second, we calculate the change 
of elastic energy due to the motion of the interface for 
fixed elastic fields in the bulk phases. After this second 
step, the coherency condition at the interface is violated. 
In the last step, we therefore have to do additional work 
to adjust the displacements appropriately. 

The first contribution is 



dW x 
dt 



dt 



(T + U (2) )dV 



y( 2 )(t) 



piHUidV -r i u lk , ik 
y<!)(t) 

(2) .(2) „r 



V( 2 'l(t) 
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We note that the kinetic energy density is continuous 
across the interface. Furthermore, by the equations of 
motion (IA5D 



dWi 
dt 



da 



(i) 

lk dV 



da 



v^'W v< 2 >(t) 



(2) 

ik dV 



d 

dxk 



dV 



(-SS- j 

y(i)(t) 



(9a 



(i) 



_9 

dx 



k v 



(2) 



W 2 )(t) 



y( 2 )(t) 



^(l) / ct^Lu^t = 0, 



A(t) 



A(t) 



where we assumed for simplicity that Ui — on all bound- 
aries apart from A(t), i. e. no external work is exerted to 
the solids. In the last step, we used the boundary condi- 
T ' ' ' — cr!fa — —a^ 2 \ 2) ] also, by definition, the 



tions (IA6I). a 



displacement rate Ui is continuous across the interface. 
The above result is quite clear since the elastodynamic 
time evolution is purely conservative. 

The second contribution arises due to the motion of the 
interface for fixed elastic fields. We extend the elastic 
state of the growing phase analytically into the newly 
acquired region. This assures that the bulk equations 
remain fulfilled in both phases even after the forward 
motion of the interface. Thus this contribution to the 
energy change rate reads 



dW 2 
dt 



A(t) 



The interface normal velocity is positive if the phase 1 lo- 
cally extends. Here, we immediately used the continuity 
of the kinetic energy density, which therefore cancels. 

After the phase transformation in this second step, the 
displacements are no longer continuous at the interface. 
Thus extra work has to be invested to remove this misfit. 
In the local coordinate system n and r (see Fig. Q} the 
strain tensor becomes 



: e rn — ^ (9rU n + 8 n U T - KU T ) . 



(A8) 
(A9) 

(A10) 



Here, k is the interface curvature, which is positive if the 
phase 1 is convex. 

At this point, a few comments concerning the conti- 
nuity of various fields across the coherent interface are 
in order. Since the displacement field has to be con- 
tinuous across the interface, also its tangential deriva- 
tives are continuous, but the normal derivatives are 
not. Consequently, the following quantities are con- 
tinuous: d T u T ,d T On the other hand, 



d n u n , d n u T , e nn , e nT are discontinuous across the inter- 
face. 

In the second step of energy calculation, we extended 
smoothly the fields into the receding domain. The inter- 
face at this new time t + At is now located at a different 
position. This leads to discontinuities of the displace- 
ments, e. g. for the normal component at the new position 
of the interface 



Au r , 



[d n u n ]^ - [cU„P>) v n At = (e$ - e™)v n At, 



where [. . denotes the evaluation of a probably dis- 
continuous expression at the previous interface position, 
taken for the phase a. Similarly, for the tangential com- 
ponent 

Au T = (2eW-[d T u n ]^ + [Ku T ]^~2eW + [d T u n ]W 
-[KU T ] {2) )v n At 
= 2(e$-eW)v n At. 

To zeroth order in At, the stresses at the new inter- 
face position are equal on both sides and identical to the 
stresses at the previous interface position. To reconnect 
the displacements, we have to apply the coherency work 
rate 



dW 3 
dt 



_f e (i) _ e ( 2 ) ) a _2fe (1) -e (2) W 



dr. 



A(t) 



Altogether, the change of the energy is given by 



(A7) aW 



dW _ d(W 1 + W 2 + W 3 ) 
~dT ~ dt 



A(t) 



CT (l) e (l) 



: a (D e (i) _ CT (D e (i) 



_ | I (T (2) e (2)_I (T (2) e (2)_ (T (2) (2) 

1 2 TT TT Q nn nn nr TIT 



dr. 



We can therefore define an appropriate chemical poten- 
tial for each phase at the coherent interface 

,/ Q ) - of _ A") _ a( a ) f ( a A CAin 

Notice that, in contrast to a free surface, the normal and 
shear contributions appear with negative sign. Then, the 
energy dissipation rate can be written as 



— = - [ (u {1) -U {2) )v dr 

A(t) 



(AI2) 



Due to the coherency condition and the requirement of 
equal mass density, the kinetic energy density does not 
appear. 
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strip (per unit length) as function of L^ 1 ': 
1 «5 2 (2 M ( 1 )+A( 1 )) 



L\ > + (L - Ly >) 2 ( 2)+A(2) 



(B3) 




FIG. 8: Motion of a planar interface in a strip geometry. 
Vertically, a constant displacement S is applied to the grips 
of the strip which consists of two solid phases with different 
elastic constants. 



This function is monotonic in Assuming that phase 
1 is harder than phase 2, 2/iW + A (1) > 2^ + A (2) , the 
energy is minimized if — 0, i.e. if the hard phase 
disappears. Notice that on the other hand, the elastic 
energy density is higher in the softer phase, < 
ojfcMfcV 2 - We get from Eq. (|B3| 



APPENDIX B: MOTION OF A PLANAR PHASE 
BOUNDARY 

The elastic contribution to the chemical potential (J5|) 
at the interface between two solid phases has the remark- 
able property that the normal and shear contributions are 
negative definite. This means that growth of the phase 
with higher elastic energy density at the expense of the 
phase with lower energy density can still reduce the total 
energy. In order to illuminate this point, we consider a 
simple example. 

Here, two strips of different solid materials are coher- 
ently connected (see Fig. [5]), and the interface can move 
due to a phase transition. We assume for simplicity 
that the process is slow and inertial effects can be ne- 
glected. We apply a fixed displacement 5 at the upper 
end of this layer structure and set the displacement at 
the lower grip to zero. The total strip width L is dis- 
tributed among the two layers, L — £,W + = const. 
We have a homogeneous strain situation in each of the 



phases (a = 1,2), i.e. u 



(a) _ 



0, u 

(a) 



(«) _ „» 



u y (y) with 



a pure linear dependence of u y on y. Strains are 



>) 

t XX 



>) 

t X y 



0, e yy 



OyUy 



const. The elongation 

of each phase is — L^e^ . In sum, they are equal 
to the prescribed total opening 5 — 5^ + = const. 
Stresses and strain are connected through Hooke's law 
for isotropic materials, thus a^x = X^e V y, o^xy = 0, 
cr$ = (2// Q ) + A (Q ))e^ ) . At the interface, the equality 
of normal stresses, a 



(i) 
yy 



(2) 

Oyy , leads to 



(2) 2/xW+AW m 

yy 2^( 2 ) + a( 2 ) yy ' 



(Bl) 



The strain can be computed in terms of the given total 
opening, 



yy 



r(i) i /r wm V^+A' 1 ' 



(B2) 



Hence, we can calculate the total elastic energy in the 



2 yy yy 



2 yy yy 



9. 



u (2) 



Here we clearly see that the negative elastic energy den- 
sity for the normal direction enters into the energy change 
rate and thus into the chemical potential; in more general 
cases one can easily verify that this is also true for shear 
contributions. 



APPENDIX C: PHASE FIELD MODELING OF 
COHERENT SOLIDS 



The aim of this section is to show that the phase field 
model presented in Section IIVI leads to the correct sharp 
interface limit even for solid-solid transformations, where 
the chemical potential differs from the usual expression 
[H[ . Since the treatment of the surface energy contribu- 
tion is well-known and enters additively into the chemical 
potential, we focus on the elastic fields here. 

For smooth interfaces, all stresses remain finite even 
in the sharp interface limit. We note that this statement 
holds even for fracture processes, where usually stresses 
can diverge at sharp corners and tips in the framework of 
the linear theory of elasticity; nevertheless, in the current 
description, the tip radius ro is always finite and there- 
fore stresses are limited to values a ~ Kr Q 1 , where 
K is a stress intensity factor. Consequently, the dis- 
placement field must be continuous in the sharp inter- 
face limit, because a finite mismatch 5u would lead to 
diverging stresses and strains e ~ a/E ~ 5u/£, and thus 
a divergent elastic energy. Then, obviously, the equation 
of motion (fT5)) leads to the usual elastic bulk equations 
and also to the continuity of normal and shear stresses 
at the interface in the limit £ — > 0. 

Next, we calculate the total elastic energy change due 
to the interface motion. We determine the energy contri- 
butions and changes in a line perpendicular through the 
interface, assuming the interface curvature to be small, 
K^Cl. To do so, we introduce a local coordinate n nor- 
mal to the interface (pointing from phase 1 into phase 2). 
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The total elastic energy per unit width is defined as 

oo 

W = J dn Qk 2 + f^j , 

— oo 

and we introduce the total elastic energy 

F ei = [ feldV, 



which is a functional of the displacement and the phase 
field. Thus 



dW 
~dT 



dn 



5F el . SF e i 86 . .. 

i — u i + ~rr^~ + P u i u i 

OUi 



56 at 



(2) 



dn 



SF el 



at 



(i) 



The first and the last term cancel each other due to the 
equation of motion f) : pure elasticity conserves en- 
ergy. Hence, the integration interval can be restricted 
to a thin region around the interface, since dissipation 
occurs only here; in fact, dt4> decays exponentially on 
the scale £. The limits of integration are inside the bulk 
phases 1 and 2, i.e. a few interface widths £ away from 
the transition point. This corresponds to the region of 
"inner equations" , as it is typically considered for rigor- 
ous sharp interface calculations of phase field models, see 
e.g. We note that upon reduction of £, the length of 
the integration interval becomes smaller proportionally. 

We furthermore assume that the interface profile moves 
without shape changes, i.e. dt = —v n d n . Then we have 



(2) 

— = -v I dn SFel 9<t> 
dt " J 56 dn 

(i) 



(CI) 



It gives explicitly 

(2) 

dW f 

— = -v n dn 



(2) 



du 2 1 d\ 2 
dnl^ k + 2dn~ eLl 



(i) 

(2) 

+ 2v n dnUe lk — + 2^—1 



(i) 



2,^2 



1(2) 
(1) 



(2) 



i , i 9e nn de nT 
+ v n dn\o nn — \-2<r nT — — 



de TT 
dn 



with the local stress a 



2u((t))e tj + \(4>)eii5ij /2. As 



we have seen above, e TT is continuous across the inter- 
face in the sharp interface limit, and therefore d n t T T is 
finite (it only has a kink at the interface for £ — * 0). 
In the sharp interface limit, the integration interval be- 
comes infinitely small, and therefore the last term in the 
integral, a TT d n e TT , does not contribute since it remains 
finite. The other terms behave differently: The stress 
components a nn and a nT are even continuous due to the 
boundary conditions. Since they vary only smoothly on 
the integration interval, they can be taken out of the inte- 
grals in the sharp interface limit. However, e nn and e nT 
are already discontinuous, and therefore, their normal 
derivatives contain delta function spikes at the interface. 
Thus, integration gives in the limit £ — > e.g. for the 
normal stress contribution 



(2) 



de nn 
1 dn 



(2) 



dn 



de r . 



dn 



■dn 



(2) 



(i) 



(i) 



Hence, we finally obtain 

dW -Vn , 
~dt~ ~ ~Y 



' nn^nn 



2a r , 



l(2) 
HI) 



■ — (u (2) ~u {1) ) 



(C2) 



with the chemical potential ((6|). 

On the other hand, we know that a solution of the 
phase field equations for an almost straight static in- 
terface is 4>{n) = (1 — tanh(n/£))/2. Here, the value 
6=1 corresponds to the phase 1. In the spirit of a rig- 
orous sharp interface analysis, where driving force terms 
behave as perturbations, this is replaced by 6{n, t) = 
(1 — tanh[(n — v n t)/(])/2 if the interface starts to move, 
here due to elastic forces. We insert the equation of mo- 
tion (f20|) into the dissipation rate (|C1[) and note that 
the double well potential and the surface energy do not 
contribute to the energy dissipation for a flat interface. 
Assuming steady state motion, we obtain 



dW 
~dt 



D 



— ) dn. 

dn, 



(C3) 



Using the above phase field profile, this gives 

oo 

dW 37V 2 , f .„ , 2 ra.o , u 2 7 

= — / (1 - tanh 2 - Ydn = — —. (C4) 

dt A£,D J V D 



Comparison with (|C2[) hence gives the sharp interface 
limit 



- — (u {1) -u (2) ) 



(C5) 



(i) 



An additional curvature contribution due to surface en- 
ergy then leads to Eq. ([7]). 
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FIG. 9: Interface velocity v/vr for the planar front scenario as 
depicted in Fig. [8] obtained from phase field simulations. The 
kinetic coefficient D is small and thus the velocity remains 
far below the Rayleigh speed. Hence, the comparison to the 
quasistatic prediction Eq. ()C5|) leads to a good agreement if 
the separation of the scales L/£ is improved. In particular, 
the interface moves such that the soft phase extends. 



We also checked this scenario numerically. For low 
kinetic coefficients D, the propagation is slow, and the 
dynamical code reproduces static elasticity. Notice that 
the sign inversion for the normal and shear terms in the 
chemical potential leads to a growth of the softer phase. 
We see the same behavior in the phase field simulations, 
and the front velocity is plotted in Fig. With increas- 
ing separation of the system size L in comparison to the 
interface width £, the interface velocity approaches the 
theoretical prediction (|C5[) . 



APPENDIX D: IMPLEMENTATION OF THE 
PHASE FIELD MODEL 

In this section, we explain in more detail the numer- 
ical discretization procedure which is designed to ob- 
tain a stable numerical algorithm for the elastic problem 
with moving boundaries. For simplicity, we use explicit 
schemes for both the phase field and the elastic equa- 
tions of motion. The dissipative phase field dynamics is 
rather robust and therefore we do not explain the pro- 
cedure here. In contrast, the elastic equations of motion 
conserve energy, and tiny numerical errors can therefore 
easily destroy the solution. We point out that energy 
conservation follows from the continuous time translation 
symmetry which is violated in any numerical discretiza- 
tion approach. Therefore, at least fluctuations in energy 
are natural, but it has to be assured that the average 
energy does not change in time. We experienced that 
naive discretization procedures can lead to long time in- 
stabilities. The generic approach which we present here is 
symmetric in time and does not suffer from this problem. 
It is not specifically related to the phase field description 
and can easily be extended to three dimensional systems 
or spatially varying mass densities. 




FIG. 10: The staggered grid: Shear modulus fi and Lame 
coefficient A are defined on the nodes (circles), the displace- 
ments Ui on the connecting lines. Thus we have three different 
lattices which are shifted by Ax/2. 



We do not discuss boundary conditions and concen- 
trate on bulk properties here. The equation of motions 
can be obtained from variational principles, as was al- 
ready shown in the preceding part of the manuscript. 
The elastodynamic evolution Eq. (fTfj|) follows from the 
action Eq. (|A4|) 



5 a. 



= 0. 



We elaborate the contributions from the kinetic and the 
potential energy separately: 



Si 



1 



-pill dV dt, 



S, 



u ■- 



1 



:<JijHij dV dt, 



and obtain for the potential part 



Su 



-5//K 2 " 

+Afi4 y ]dVdt. 



We use a staggered grid, i.e. the mass density and the 
elastic constants are defined on the grid points, displace- 
ments between them (see Fig. [Tt)|) [24j . In our case, the 
spatial (and temporal) values of the elastic coefficients 
/i, A are related to the phase field. Similar to the deriva- 
tion above, we keep the phase field fixed (and thus the 
elastic coefficients) during the variation with respect to 

the elastic displacements. We use the notation u^\i,j), 
where i, j are the spatial and n is the time index; in the 
phase field formulation, no explicit distinction between 
the different phases has to be made, and therefore the 
upper index cannot be confused with previous notations. 
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We assume the grid spacing Ax to be the same in both 
spatial directions. 

The central idea for derivation of the discrete equations 
of motion is the discretization of the action (obeying sym- 
metry in space and time) and to perform discrete vari- 
ations with respect to each degree of freedom u x n \i,j) 
and u y n \i,j). We study the potential contribution to S 
first: 



5(7 



-\(Ax) 2 At 



^xx 

+ 6y y ) + 2\e xx e yy 



n i,j 



on grid points 



+ 



^^xy 



in square center 



We express the first part on the grid points, and therefore 
replace the elastic coefficients as follows: 

M^M(^i), A->A(i,j). 
Strains also have to be evaluated on the nodal points: 



e - c (") (i i) - u[ - ] ^3)-u x n, {i-l,3) 



An) 

yy 



u y ihJ) ~ Uy (hJ ~ 1) 
Ax 



The second part is expressed in the center of the squares, 
i. e.: 

p -» a*(* + 1/2,J + 1/2) 

ex, - e<#(i + l/2,j + l/2) 

4")(i,i + l)-4")(i,i)+4")(i + l,i) 

j)]/(2Aa;). 

We illustrate the discrete variation with respect to 

(n) / . - \ 



OSu 



dui n \i,j) 



-AtAx{[2p(i,j) + \(i,j)] e x n J(i,j) 



-[2p(i + l,j) + \(i + l,j)]e x n J(i + l,j) 



+\(i,j)e y n y >(i,j)-\(i + l,j)e y n y >(i + l,j) 
-2/x(t + 1/2, j + l/2)eW(i + 1/2, j + 1/2) 
+2»(i + 1/2, j - l/2)e x n y \i + 1/2,3 - 1/2)}. 

For the kinetic contribution, we proceed in a similar way. 
Here, the terms are defined between the lattice points: 



5 T ^^(Ax) 2 A^^( pii 2 x + P il 



Discretization of the first term defines the displacement 
rate v x n+1 ^ 2 \i, j) at intermediate timesteps 



At 



(»)/ 



and similarly for the second term 



.,„+! 2),; ._ U { y +1 \i,i) -Uy n \i,j) 



Uy^vW*>(i,j) := 



At 



Variation of the kinetic contribution to the discrete action 
therefore gives 



dui n \i,j) 



-(Ax) 2 pAt 



,x n+1) ^3)-2u x n \ l ,j) + ux n - 1 \ l ,j) 
(At) 2 



Notice that this expression is invariant against time in- 
version. Vanishing total variation of S — Sjj + St with 
respect to u x n \i,j) leads to the desired explicit evolution 
equation. 

The same procedure has to be performed for u y . 

We performed various tests to check the code, among 
them the verification of the sound speeds. The theo- 
retical expressions for the dilatational and shear wave 
speed, Vd = [(A + 2p)/p] 1 / 2 and v s = (p/p) 1 ^ 2 were ob- 
tained with high accuracy. Also, we checked the trans- 
mission and reflection coefficients of both wave types at 
stationary interfaces. Here, we froze the dynamics of the 
phase fields and let shock waves hit the straight interface 
between different solid phases. The impedance of each 
phase is defined as = pv^ with the relevant sound 
speed for the considered wave type in each phase a. 
The reflection coefficient R is defined as the ratio of the 
amplitudes of the reflected and the incoming wave and is 
given by 



R = 



(Dl) 



and the transmission coefficient is similarly T = 1 + R. 
Both values are reproduced by the numerical simulations 
for a phase boundary between two solids. 
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